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Abstract 

A new approximate Riemann solver for the equations of magnetohydro dynamics 
(MHD) with an isothermal equation of state is presented. The proposed method of 
solution draws on the recent work of Miyoshi and Kusano, in the context of adia- 
batic MHD, where an approximate solution to the Riemann problem is sought in 
terms of an average constant velocity and total pressure across the Riemann fan. 
This allows the formation of four intermediate states enclosed by two outermost 
fast discontinuities and separated by two rotational waves and an entropy mode. 
In the present work, a corresponding derivation for the isothermal MHD equations 
is presented. It is found that the absence of the entropy mode leads to a different 
formulation which is based on a three-state representation rather than four. Numer- 
ical tests in one and two dimensions demonstrates that the new solver is robust and 
comparable in accuracy to the more expensive linearized solver of Roe, although 
considerably faster. 
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1 Introduction 

The numerical solution of the magnetohydrodynamics (MHD) equations has received 
an increasing amount of attention over the last few decades motivated by the growing 
interests in modeling highly nonlinear flows. Numerical algorithms based on upwind dif- 
ferencing techniques have established stable and robust frameworks for their ability to 
capture the dynamics of strong discontinuities. This success owes to the proper compu- 
tation of numerical fluxes at cell boundaries, based on the solution of one-dimensional 
Riemann problems between adjacent discontinuous states. The accuracy and strength of 
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a Riemann solver is intrinsically tied to the degree of approximation used in capturing 
the correct wave pattern solution. 

In the case of adiabatic flows, an arbitrary initial discontinuity evolves into a pattern 
of constant states separated by seven waves, two pairs of magneto-sonic waves (fast and 
slow), two Alfven discontinuities and one contact or entropy mode. With the exception 
of the entropy mode, a wave separating two adjacent states may be either a shock or a 
rarefaction wave. In the limit of efficient radiative cooling processes or highly conducting 
plasma, on the other hand, the adiabatic assumption becomes inadequate and the ap- 
proximation of an isothermal flow is better suited to describe the flow. In this case, the 
solution of the Riemann problem is similar to the adiabatic case, with the exception of 
the contact mode which is absent. 

Although analytical solutions can be found with a high degree of accuracy using exact 
Riemann solvers, the resulting numerical codes are time-consuming and this has pushed 
researcher's attention toward other more efficient strategies of solution. From this perspec- 
tive, a flourishing number of approaches has been developed in the context of adiabatic 
MHD, see for example [4a-b|TTTll3ll 1 a-b] and references therein. Similar achievements have 
been obtained for the isothermal MHD equations, see [la-b|[TU] . 

Recently, Miyoshi and Kusano ([9], MK henceforth) proposed a multi-state Harten-Lax- 
van Leer (HLL, [5]) Riemann solver for adiabatic MHD. The proposed "HLLD" solver 
relies on the approximate assumption of constant velocity and total pressure over the Rie- 
mann fan. This naturally leads to a four-state representation of the solution where fast, 
rotational and entropy waves are allowed to form. The resulting scheme was shown to be 
robust and accurate. The purpose of the present work is to extend the approach devel- 
oped by MK to the equations of isothermal magnetohydrodynamics (IMHD henceforth). 
Although adiabatic codes with specific heat ratios close to one may closely emulate the 
isothermal behavior, it is nevertheless advisable to develop codes specifically designed for 
IMHD, because of the greater ease of implementation and efficiency over adiabatic ones. 

The paper is structured as follows: in $2] the equations of IMHD are given; in $3] the new 
Riemann solver is derived. Numerical tests are presented in £JH and conclusions are drawn 
in 83 



2 Governing Equations 



The ideal MHD equations describing an electrically conducting perfect fluid are expressed 
by conservation of mass, 

^ + V-(pv) = 0, (1) 

momentum, 

— + V • (mv - BB + Ip T ) = , (2) 

and magnetic flux, 

f)T2 

— -Vx(vxB) = 0. (3) 
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The system of equations is complemented by the additional solenoidal constraint V-B = 
and by an equation of state relating pressure, internal energy and density and Here p, v 
and m = pv are used to denote, respectively, density, velocity and momentum density 
The total pressure, px, includes thermal (p) and magnetic (|B| 2 /2) contributions (a factor 
of v^47r has been reabsorbed in the definition of B). No energy equation is present since 
I will consider the isothermal limit for which one has p = a 2 p, where a is the (constant) 
isothermal speed of sound. With this assumption the total pressure can be written as: 

IBI 2 

PT = a 2 p+ 1 -^-. (4) 



Finite volume numerical schemes aimed to solved ([I])-© rely on the solution of one- 
dimensional Riemann problems between discontinuous left and right states at zone edges. 
Hence without loss of generality, I will focus, in what follows, on the one- dimensional 
hyperbolic conservation law 



d\J dF _ 

dt dx 



(5) 



where the vector of conservative variables U and the corresponding fluxes F are given by 
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(6) 



where u, v and w are the three components of velocity in the x, y and z directions, 
respectively. The solution of Eq. (jSj) may be achieved through the standard two-point 
finite difference scheme 



U 



n+l 



At 
Ax 



1 2 



(7) 



where F i± i is the numerical flux function, properly computed by solving a Riemann 
problem between Uj and Uj±i. 

As already mentioned, the solution to the Riemann problem in IMHD consists of a six wave 
pattern, with two pairs of magneto-sonic waves (fast and slow) and two Alfven (rotational) 
waves. These families of waves propagate information with characteristic signal velocities 
given by the eigenvalues of the Jacobian <9F(U)/<9U: 



A 



.f± 



U± Cf , A; 



± = u ± c s , A a± = u ± c a , 



(8) 
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Fig. 1. Schematic representation of the approximated wave structure used to describe the Rie- 
mann fan. The initial states Ul and Ur are connected to each other through a set of four waves 
(Sl, S l S r and Sr) across which some or all of the fluid variables can be discontinuous. The 
four discontinuities isolate three intermediate constant states J7£, U* and U R which correspond 
to the approximate solution of the Riemann problem. 

where the fast, Alfven and slow velocity are denned as in adiabatic MHD: 
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1/2 



(9) 



with b = ~B/y/p. 
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(10) 

(11) 



The contact or entropy mode is absent. Fast and slow waves can be either shocks, where 
flow quantities experience a discontinuous jump, or rarefaction waves, characterized by a 
smooth transition of the state variables. Across the rotational waves, similarly to the adi- 
abatic case, density and normal component of velocity are continuous, whereas tangential 
components are not. In virtue of the solenoidal constraint, the normal component of the 
magnetic field is constant everywhere and should be regarded as a given parameter. 



3 The HLLD Approximate Riemann Solver 



In what follows, I will assume that the solution to the isothermal MHD Riemann problem 
can be represented in terms of three constant states, XJ* L , U* and XJ* R , separated by four 
waves: two outermost fast magneto-sonic disturbances with speeds Sl and Sr (Sl < Sr) 
and two innermost rotational discontinuities with velocities S L and S R , see Fig. [TJ Slow 
modes cannot be generated inside the solution. Since no entropy mode is present and 
rotational discontinuities carry jumps in the tangential components of vector fields only, 
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it follows that one can regard density, normal velocity and total pressure as constants over 
the fan. This assumption is essentially the same one suggested by MK but it leads, in the 
present context, to a representation in terms of three states rather than four. Quantities 
ahead and behind each wave must be related by the jump conditions 

S a (U* a -U a ) = F* a -F a , (12) 

across the outermost fast waves, and 

s; (u; - u;) = f; - f; , (is) 

across the innermost Alfven waves. Here a = L or a = R is used for the left or right state, 
respectively. Note that U* and F* are both unknowns of the problem and F* should not 
be confused with F(U*), although one always has F a = F(U Q ). 

The solution to the problem must satisfy the consistency condition which, in terms of U*, 
can be obtained by direct summation of Eq. (fl2"l) and ffTB"]) across all waves: 

(si - s L ) vi + (s R - sp u; + (s R - s* R ) u| = uhll ^ 

Sr — Si 

where 

Uhll _ SrUr ~ S l Ul - F R + F L ^ 
Sr — Si 

is the HLL single average state. Likewise, after dividing Eq. (fl2"|) and (|T3|) by S a and S* 
(respectively) and adding the resulting equations, the following condition for the fluxes 
F* is found: 

1 1 \ Tp* I / 1 1 \ TP* I I 1 1 I TP* 

- - L ) F ^ + \m - si) F c + [si - s% f r , u 



j i_ 

Sr Sl 



where 



phll 



SrFl — SlFr + SrSl(Ur — Ul) 



F mi , (16) 



(17) 



Sr. — Sl 

is the HLL single average flux. Equation (1171) represents the average flux function inside 
the Riemann fan and can be considered as the "dual" relation of (|15)) . Indeed, Eq. (fTTj) 
could have been obtained by letting XJ a — >• F a , F a — >• XJ a and S a — > l/S a on the right 
hand side of Eq. ffT5l) . 

In this approximation, density, normal velocity and total pressure do not experience jumps 
across the inner waves and are constant throughout the Riemann fan, i.e., p* L = p* = p* R = 
p*, u* L = u* c = u* R = u*, Pt l = Vt c = Pt r = Pt- Tangential components of magnetic field 
and velocity, on the other hand, may be discontinuous across any wave and their jumps 
must be computed using (|I2"1) or (fT3"j) . 

At first sight, one might be tempted to proceed as in the adiabatic case, that is, by writing 
U* and F* in terms of the seven unknowns p*, u*, v*, w*, B*, B* and p?, and then using 
Eq. ( JT2|) and (fT3"j) . together with the continuity of p*, u* and p?, to determine the jumps 
accordingly. However, this approach is mathematically incorrect since the resulting system 



5 



of equations is overdetermined. This can be better understood by explicitly writing the 
first and second component of Eq. i\V2\i : 



S a (p* - Pa) = P*U* ~ PaUa , (18) 

S a (p*u* - p a u a ) =p*(u*) 2 + p* T - p a u a -p Ta , (19) 

where I have used the fact the p*, u* and p* T are constants in each of the three states. 
Eqns ( I18p and f[T9"j) provide four equations (two for a = L and two for a = R) in the 
three unknowns p*, u* and p? and, as such, they cannot have a solution. Furthermore, 
the situation does not improve when the equations for the tangential components are 
included, since more unknowns are brought in. 

The problem can be disentangled by writing density and momentum jump conditions in 
terms of four unknowns rather than three. This is justified by the fact that, according to 
the chosen representation, density and momentum components of U* and F* are contin- 
uous across the Riemann fan and, in line with the consistency condition, they should be 
represented by their HLL averages: 

* * * * till /nn\ 

P = P L = Pc= PR = P , (20) 



* * * # hll /ni \ 

m x = m x L = m x c = m x R = m x , (21) 



17* TP* JP* TP* I?hll fOO\ 

*\p] = *\p]l - h W - h \p] R ~ *\p\ ' (22) 



TP* — TP* TP* Tp* Z?U1 (OQ\ 

r [m x ] = r [m x ]L - r [rn x }c ~ 1 [m x ]R ~ r [m x ] ) K A6 ) 

where p hU , m hU , and j are given by the first two components of Eq. (Tl5|) and Eq. 

(El)- 

This set of relations obviously satisfies the jump conditions across all the four waves. 
The absence of the entropy mode, however, leaves the velocity in the Riemann fan, u*, 
unspecified. This is not the case in the adiabatic case, where one has m* = Fj*i = p*u* 
and from the analogous consistency condition the unique choice p hll u* = m^ 11 . Still, in 
the isothermal case, one has the freedom to define u* = F^/p hYl or u* = m^ u /p hU . The 
two choices are not equivalent and one can show that only the first one has the correct 
physical interpretation an "advective" velocity. This statement becomes noticeably true 
in the limit of vanishing magnetic fields, where transverse velocities are passively advected 
and thus carry zero jump across the outermost sound waves. Indeed, by explicitly writing 
the jump conditions for the transverse momenta, one can see that the conditions v * = v a 
and w* a = w a can be fulfilled only if p*u* = Ffa . 
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Thus one has, for variables and fluxes in the star region, the following representation: 
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which gives 8 variables per state and thus a total of 24 unknowns. The 4 continuity 
relations (T20T) — (T23T) together with the 12 jump conditions (JT2"j) across the outer waves can 
be used to completely specify U* (for a = L and a = R) in terms of Ul and Ur. Note 
that the total pressure does not explicitly enter in the definition of the momentum flux. 
The remaining 4 unknowns in the intermediate region (v*, w*, B* c and B*J may be found 
using (fl3|) . However, one can see that 



(s* a - u*) P *v: + b x b; c = (s* a - «>x + b x b^ , (25) 

(S* a - u*) P *w* c + B x B* Zc = (S* a - u*)p*w* a + B x B* Za , (26) 

(S* a - u*)B* yc + B x vl = (S* - u*)B* ya + B x v* a , (27) 

(S* a - u*)B* Zc + B x w* c = (S* a - u*)B* Za + B x w* a , (28) 

cannot be simultaneously satisfied for a = L and a = R for an arbitrary definition of 5*, 
unless the equations are linearly dependent, in which case one finds 

Sl = u*- l -^4, S* R = u*+ l -^4. (29) 



This choice uniquely specifies the speed of the rotational discontinuity in terms of the 
average Alfven velocity \B x \/y/ff. Thus only 4 (out of 8) equations can be regarded 
independent and one has at disposal a total of 24 independent equations, consistently 
with a well-posed problem. 



291) and the assumption u* = F<^/ p m yields 



Proceeding in the derivation, direct substitution of (12411 into Eqns (112]) together with (1201) . 

Tlhll 

[p] 



P*K = P*Va - B x B Va — a — — , (30) 

u* u 
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b: 



which clearly shows that {v* a ,w* ai 



B Za Pa {S a -u a f-Bl 
P* (S a — S1)(S a — S R ) 



(33) 



(v a , w a ) when B 



0. As in the adiabatic case, the 



consistency condition should be used to find v*, w*, B* and B* : 



P v, 



* * 

P W c 



( P *wi + p*w* R ) | x(b; r -b* zl 



Vc 2 2X 

B * = (B* ZL + B; R ) (p*w* R - p*w* L ) 



2X 



(34) 
(35) 
(36) 
(37) 



where X = y/ffsign(B x ). The inter-cell numerical flux F required in Eq. (J5j) can now be 
computed by sampling the solution on the x/t = axis: 



F L for S L > 

F L + S L (V* L - U L ) for S L < < S* L 

K ^ S* L <0<S R , (38) 

F fi + S R (V* R - U fl ) for 5^ < < S R 

F r for S R <0 



where F* is computed from ([24"]) using (|34j) through (1371) . 

To conclude, one needs an estimate for the upper and lower signal velocities Sl and S R . 
I will adopt the simple estimate given by [5], for which one has 



S L = min (A/_(Ux), A/_(Ur)) , = min (A /+ (U L ), A /+ (U R )) 



(39) 



with A/±(U) given by the first of ([8]). Other estimates have been proposed, see for example 
0. " ' 

The degenerate cases are handled as follows: 

• for zero tangential components of the magnetic field and B 2 X > a 2 p a , a degeneracy 
occurs where 5* — ► S a and the denominator in Eqns. (I30l) - fl33l vanishes. In this case, 
one proceeds as in the adiabatic case by imposing transverse components of velocity 
and magnetic field to have zero jump; 

• the case B x — >• does not poses any serious difficulty and only the states and XJ* R 
have to be evaluated in order to compute the fluxes. One should realize, however, that 
the final representation of the solution differs from the adiabatic counterpart. In the 
full solution, in fact, Alfven and slow modes all degenerate into a single intermediate 
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tangential discontinuity with speed u*. In the adiabatic case, this degeneracy can still 
be well described by the presence of the middle wave, across which only normal velocity 
and total pressure are continuous. Although this property is recovered by the adiabatic 
HLLD solver, the same does not occur in our formulation because of the initial assump- 
tion of constant density through the Riemann fan. The unresolved density jump, in 
fact, can be thought of as the limiting case of two slow shocks merging into a single 
wave. Since slow shocks are not considered in the solution, density will still be given by 
the single HLL-averaged state. Tangential components of velocity and magnetic field, 
on the other hand, can yet be discontinuous across the middle wave. Thus, in this limit, 
the solver does not reduce to the single state HLL solver. 

This completes the description of the proposed isothermal HLLD Riemann solver. The 
suggested formulation satisfies the consistency conditions by construction and defines a 
well-posed problem, since the number of unknowns is perfectly balanced by the number of 
available equations: 12 equations for the outermost waves, 8 continuity conditions (for p*, 
u*, Fr*,, KJ^ j across the inner waves) and the 4 equations for the tangential components. 
The positivity of the solver is trivial since no energy equation is present and the density 
is given by the HLL-averaged state (!20|) . 



4 Numerical Tests 



The performance of the new Riemann solver is now investigated through a series of se- 
lected numerical tests in one and two dimensions. One dimensional shock tubes (in §4.11) 
consist in a decaying discontinuity initially separating two constant states. The outcoming 
wave pattern comprises several waves and an analytical or reference solution can usually 
be obtained with a high degree of accuracy. For this reason, one dimensional tests are 
commonly used to check the ability of the code in reproducing the correct structure. 

The decay of standing Alfven waves in two dimensions is consider in §4.21 The purpose 
of this test is to quantify the amount of intrinsic numerical dissipation inherent to a 
numerical scheme. A two dimensional tests involving the propagation of a blast wave in a 
strongly magnetized environment is investigated in §4.31 and finally, in §4.41 the isothermal 
Orszag-Tang vortex system and its transition to supersonic turbulence is discussed. 

Eq. ([7]) is used for the spatially and temporally first order scheme. Extension to second 
order is achieved with the 2 nd order fully unsplit Runge-Kutta Total Variation Diminishing 
(TVD) of [7] and the harmonic slope limiter of [H]. For the multidimensional tests, the 
magnetic field is evolved using the constrained transport method of [2], where staggered 
magnetic fields are updated, e.g. during the predictor step, according to 
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Table 1 

Initial conditions for the one-dimensional shock tube problems presented in the text. Here 
(H x , H y , H z ) = \/ r 4^r(B x , B y , B z ). The final integration time (t s ) is given in the last column. 

The electric field Vt z i+ 1 . +1 is computed by a spatial average the upwind fluxes available 
at zone interfaces during the Riemann solver step: 

p[By\ p[By\ _|_ p[Bx] _|_ p[B x } 

4 ' ^ > 

where F^ By ^ is the B y component of the flux ( 1381) available during the sweep along the 
x-direction. A similar procedure applies to F^ Ba \ 



4-1 One Dimensional Shock Tubes 



An interface separating two constant states is placed in the middle of the domain [0, 1] at 
x = 0.5. Unless otherwise stated, the isothermal speed of sound a is set to unity. States 
to the left and to the right of the discontinuity together with the final integration time 
are given in Table [H The reader may also refer to [10] for the first and second problems 
and to |lb] for the third one. For the sake of comparison, the linearized Riemann solver 
of Roe (properly adapted to the isothermal case from [3]) and the simple single-state 
HLL solver of [5] will also be used in the computations. In order to better highlight the 
differences between the selected Riemann solvers, integrations are carried with a spatially 
and temporally first order-accurate scheme on 400 uniform zones. The Courant number 
is set to 0.8 in all calculations. 

Profiles of density, velocities and magnetic fields for the first shock tube are shown in Fig. 
([2]) . Solid, dashed and dotted lines refer to computations carried out with the HLLD, 
Roe and HLL solvers, respectively. The decay of the initial discontinuity results in the 
formation, from left to right, of a fast and a slow rarefaction waves followed by a slow 
and a fast shocks. No rotational discontinuity is present in the solution. The leftmost 
rarefaction fan and the rightmost fast shock are equally resolved by all three solvers. 
The resolution of the slow rarefaction is essentially the same for the Roe and the HLLD 
schemes, but is smeared out over more computational zones for the HLL method. At 
the slow shock, differences become less evident although the solution obtained with the 
HLLD solver seems to spread the discontinuity over a slightly higher number of zones 
with respect to the Roe scheme. Although the HLLD solver should, in principle, offer no 
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Fig. 2. Computed profiles for the first shock tube problem at time t = 0.1. The top panels 
show, from left to right, density (p), and transverse components of magnetic field {B y and B z , 
respectively). Shown on the bottom panels are the three components of velocity. Solid, dashed 
and dotted lines refer to computations carried with the HLLD, Roe and HLL Riemann solvers, 
respectively. 
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Fig. 3. Computed profiles for the second shock tube problem at time t = 0.2. The top panels show, 
from left to right, density (/?), y and z components of magnetic field (B y and B z , respectively). 
Shown on the bottom panels are the three components of velocity. Solid, dashed and dotted 
lines refer to computations carried with the HLLD, Roe and HLL Riemann solvers, respectively. 

gain over the simpler HLL method in absence of rotational waves, the wave patterns is 
still reproduced more accurately than the HLL scheme. 

The ability to resolve Alfven waves is better observed in the second shock tube (test 2) 
with initial conditions given in Table [TJ In this case one has the formation of three pairs of 
waves which, from the exterior to the interior of the domain, can be identified with a pair 
of fast magnetosonic shocks, a pair of rotational discontinuities and a pair of slow shock 
waves, see Fig. [31 One can note the poorer performance of the HLL solver in resolving 
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Fig. 4. CPU time (in units of the fastest run) for the HLLD (solid lines) and Roe schemes 
(dashed lines). The adiabatic cases (with T = 1.01) are overplotted using boxes. 

the rotational discontinuities which are smeared out on several zones. Likewise, the peak 
in the y component of magnetic field between the Alfven wave and the slow shock is 
considerably underestimated (> 4%). On the contrary, the Roe and HLLD schemes give 
comparable performances in terms of reduced numerical diffusion on essentially all waves, 
including the two slow shocks. The plots shown in Fig [3] show, in fact, almost complete 
overlap. In terms of efficiency, however, the new solver offers ease of implementation over 
the Roe solver which requires the characteristic decomposition of the Jacobian matrix. 
This is better quantified in Fig. HI where the two schemes are compared in terms of CPU 
execution times as functions of the mesh size. At the resolutions employed (200 to 6400), 
the isothermal HLLD is faster than the Roe solver by a factor ~ 35 4- 40%. The same 
performance is found if one compares the efficiencies of the two solvers by running the 
same test with an adiabatic equation of state (by setting F = 1.01). Furthermore, direct 
inspection of Fig. H] reveals that at lower resolution (i.e. < 600) the CPU time ratio 
adiabatic/isothermal is ~ 1.2, whereas it increases up to ~ 1.35 at higher resolutions. 
This justifies the need for codes specifically designed for IMHD, as already anticipated in 

m 

The collision of oppositely moving high Mach number flows is examined in the third shock 
tube, see also [lbj . This problem also illustrates the performances of the new solver in the 
limit of vanishing normal component of the magnetic field, B x = 0, as discussed at the 
end of $3j The wave pattern is bounded by two outermost fast shocks traveling in opposite 
directions and a stationary tangential discontinuity (in the middle) carrying no jump in 
total pressure and normal velocity, see Fig [51 The two shocks propagate at the same speed 
which can be found analytically by solving the jump conditions (|T2|) with respect to S a , the 
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Fig. 5. Computed profiles for the third shock tube problem at time t = 0.25. The top panels show, 
from left to right, density (p), y and z components of magnetic field (B y and B z , respectively). 
Shown on the bottom panels are the three components of velocity. Solid, dashed and dotted lines 
refer, respectively, to computations carried with the HLLD, Roe and HLL Riemann solvers. 



upstream density and magnetic field (the upstream velocity is zero). This yields the value 
\S a \ ~ 1.6529 and therefore predicted shock positions at ~ 0.5 ±0.4132 when t = 0.25. At 
the resolution of 400 zones, the HLL solver performs noticeably worse than the other two 
schemes, resulting in slower shock propagation speeds, hence the underestimated jumps. 
On the contrary, the HLLD and Roe solvers perform very similarly across all three waves 
and are able to reproduce the expected shock position with relative errors less than 10~ 2 . 
The spurious density spike at the center of the grid is a numerical artifact generated by 
the solver, which is trying to compensate the reduced magnetic pressure (caused by the 
spreading of B y and B z over a region of finite width across the discontinuity) by increasing 
the density in such a way that the total pressure is maintained constant. This effect is 
manifestly more pronounced for the HLL scheme yielding a peak value of ~ 3.1 and 
become less noticeable for the Roe scheme (~ 0.73). The HLLD solver does not exhibit 
this feature and yields a flat density profile, concordantly with the exact solution. 

In order to quantify the errors and check the convergence properties of the scheme, com- 
putations have been repeated at different grid sizes N x = 50, 100, 200, 400, 800, 1600 and 
the Li-norm errors have been measured against a reference solution. The reference so- 
lution was obtained with the second order version of the scheme using an adaptive grid 
with 3200 coarse zones and 4 levels of refinement, therefore corresponding to an effective 
resolution of 512, 000 uniform zones. Results shown in Fig. [6] confirm that the new solver 
reproduces the correct solution with considerably smaller errors than the simple HLL 
scheme. In the second test problem, the Roe and HLLD schemes perform nearly the same 
with density and transverse magnetic field errors being almost halved (< 4%, at the high- 
est resolution) with respect to the HLL method. In the third problem, the HLLD solver 
is able to capture the solution with even higher accuracy than the Roe method indeed 
yielding, at the highest resolution, errors of ~ 0.03% and ~ 0.07% in density and magnetic 
field (respectively) compared with 0.1% and 0.13% obtained with the Roe scheme. The 
unusually higher errors (> 10%) observed in the solutions obtained with the HLL scheme 
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Fig. 6. Percent LI norm errors (computed as £Li(q) = lO 2 /-/^ Yll=i x \Q ex ( x i) ~ l( x i)\ where q ex 
and q are the exact and numerical solutions) as function of the resolution, for the first (left 
panels), second (middle panels) and third (right panels) shock tube problems. Top and bottom 
rows give the error in density and transverse magnetic field B y , respectively. 

are accounted for by the higher density overshoot occurring at the center of the grid. 



4-2 Decay of Standing Alfven Waves 



The amount of numerical viscosity and/or resistivity intrinsic to computations carried 
on discrete grids may be estimated with a simple numerical test involving the decay of 
standing Alfven waves, [T21T0"] . In a fluid with non-zero dynamic shear viscosity \x and 
resistivity rj, waves are characterized by a decay rate given by 

where po is the background density and k is the wave number vector. The inevitable ap- 
proximations introduced by any numerical scheme lead to an effective intrinsic numerical 
viscosity (and resistivity) that can be also accounted for by the particular Riemann solver 
being used. Thus, a measure of Ta can be directly related to the "effective" \i and r\ 
inherent to the scheme. 

The initial condition consists of a uniform background medium with constant density 
p = 1 and magnetic field B = x. The x and y components of velocity are set to zero 
everywhere and the sound speed is a = 1. Standing Alfven waves are set up along the 
main diagonal by prescribing 

v z = £Ca sin (k • x) (44) 

where ca = l/\/2, k = 2n(x + y) and e = 0.1. Computations are carried over the unit 
square [0, l] 2 until t — 10 at different resolutions (16 2 , 32 2 , 64 2 , 128 2 zones) and for the 
three Riemann solvers. Periodic boundary conditions are imposed on all sides. Following 
[TO] , the decay rate corresponding to a given resolution (for < t < 10) has been derived 
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Fig. 7. Decay rate as function of resolution for the standing Alfven wave test. Left and right 
panels give the computed decay rate on a log-log scale for the first and the second order scheme, 
respectively. At a given resolution, the decay rate is measured by first fitting the peaks <5-B™ ax (t) 



of SB z (t) = JJJ B%(t)dxdy with respect to time, and then by approximating the corresponding 
with dlog (5i?™ ax (i)) jdt computed as average slope over the time integration interval. Note 
that the HLLD and Roe solvers yields essentially the same decay rate and the two lines cannot 
be distinguished. 

by fitting the peaks of the r.m.s. of the z component of the magnetic field, 5B z (t), with 
respect to time. Fig. shows the normalized decay rates for the first and second order 
schemes at different resolutions and for the selected Riemann solvers. The amount of 
numerical viscosity intrinsic to the HLLD and Roe scheme turns out to be virtually the 
same, whereas it is a factor of 3 (for the 1 st order) and ~ 2.6 (for the 2 nd order) higher 
for the HLL solver. The slope of the curves is the same for all of the methods, therefore 
confirming the selected integration order of accuracy (oc N~ Y or oc N~ 2 for the 1 st or 2 nd 
order scheme, respectively). 

4-3 Cylindrical Blast Wave 

A uniform medium with constant density p = 1 and magnetic field B = 10/V^bri fills the 
square domain [— |, |] 2 . An over-pressurized region (p — 100) is delimited by a small circle 
of radius 0.05 centered around the origin. The speed of sound is set to unity everywhere 
and free-flow boundary conditions are imposed on every side of the domain. 

Cylindrical explosions in Cartesian coordinates provide a useful benchmark particularly 
convenient in checking the robustness of the code and the response of the algorithm to 
different kinds of degeneracies. Fig. [8] shows the results obtained on 40 1 2 uniform zones 
at t = 0.09. The ratio of thermal and magnetic pressure is ~ 0.251 thus making the 
explosion anisotropic. The outer region is delimited by an almost cylindrical fast forward 
shock reaching its maximum strength on the y axis (where the magnetic field is tangential) 
and progressively weakening as y — > 0. In the central region, density and magnetic fields 
are gradually depleted by a cylindrical rarefaction wave enclosed by an eye-shaped reverse 
fast shock bounding the inner edge of a higher density ring. The outer edge of this ring 
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Fig. 8. Colored contours for the blast wave problem at t = 0.09. Density logarithm (upper 
left panel), magnetic pressure logarithm (upper right), specific kinetic energy (bottom left) and 
magnetic field lines (bottom right) are shown. 

marks another discontinuity which degenerates into a tangential one on the y axis and 
into a pure gas dynamical shock at the equator (y = 0) where B y = 0. The HLLD solver 
provides adequate and sharp resolution of all the discontinuities and the result favorably 
compares to the integrations carried out with the other Riemann solvers (not shown here) 
and with the results of llbl. 



The computational execution times for the HLL, HLLD and Roe solvers stay in the ratio 
1 : 1.21 : 2.01, thus showing that the new method is approximately 20% more expensive 
whereas the Roe solver results in computations which are almost twice as slow. 



4-4 Isothermal Orszag Tang Vortex 



The compressible Orszag- Tang vortex system describes a doubly periodic fluid configura- 
tion undergoing supersonic MHD turbulence in two dimensions. Although an analytical 
solution is not known, its simple and reproducible set of initial conditions has made it a 
widespread benchmark for inter-scheme comparison, see for example |13j . An isothermal 
version of the problem was previously considered by [lb] . In his initial conditions, however, 
the numerical value of the sound speed was not specified. For this reason, I will consider 
a somewhat different initial condition with velocity and magnetic fields prescribed as 

v = a (— sin yit + sin xy) , (45) 
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Fig. 9. Evolution of the Orszag-Tang vortex at t = 1.5 (top) and t = 3 (bottom). Density and 
magnetic pressure are shown on the left and on the right, respectively. Note the formation of a 
magnetic island in the center at t = 3. 

B = B (- sin y± + sin(2x)y) . (46) 

The value of the speed of sound is a = 2 and the magnetic field amplitude is Bq = aJ3/5. 
The density is initialized to one all over. This configuration yields a ratio of thermal 
to magnetic pressure of 10/3, as in the adiabatic case. The domain is the box [0,27r] 2 
with periodic boundary conditions imposed on all sides. Fig. O shows the density and 
magnetic pressure distributions at t = 1.5 and t = 3. The initial vorticity distribution 
spins the fluid clockwise and density perturbations steepen into shocks around t ~ 0.7. 
Because of the imposed double periodicity, shocks re-emerge on the opposite sides of 
the domain and broad regions of compression are formed. The dynamics is regulated by 
multiple interactions of shock waves, which lead to the formation of small scale vortices and 
density fluctuations. The shocks continuously affect the spatial distributions of density and 
magnetic field. Magnetic energy is dissipated through current sheets convoyed by shocks 
at the center of the domain where reconnection takes place around t ~ 1.5, through the 
"Y" point. By t = 3 the kinetic energy has been reduced by more than 50% and shocks 
have been weakened by the repeated expansions, leaving small scale density fluctuations. 

Note that a magnetic island featuring a high density spot forms at the center of the 
domain by t — 3. This corresponds to a change in the magnetic field topology where an 
"O" point forms in the central region. It should be emphasized that the same feature is 
observed when the Roe solver is employed during the integration and it is also distinctly 
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recognizable in the results of [Tbj. However, it requires at least twice the resolution to be 
seen if the fluxes are computed with the HLL solver. Yet, the magnetic island does not 
form at sufficiently lower resolution for any solver. Although its origin has to do with the 
complex phenomenology of MHD reconnection which lies outside the scope of this paper, 
it is speculated that its presence may be numerically induced by a decreased effective 
resistivity across the central current sheet. Future studies should specifically investigate 
on this issue. 

The total execution times for the three solvers (HLL, HLLD and Roe) scale as 1 : 1.25 : 
2.06, thus confirming the results obtained for the previous test problem. 



5 Conclusions 

A new Riemann solver for the equations of magnetohydrodynamics with an isothermal 
equation of state has been derived. The solver leans on a multi-state Harten-Lav-van Leer 
representation of the solution inside the Riemann fan, where only fast and rotational 
discontinuities are considered. Under this approximation, density, normal momentum and 
the corresponding flux components inside the solution are given by their HLL averages. 
This is equivalent to the assumption that normal velocity, density and total pressure may 
be regarded constant across the solution. It has been shown that this consistently leads 
to a three-state representation where only tangential vectors may be discontinuous across 
the inner rotational waves. The solver is simple to implement and does not require a 
characteristic decomposition of the Jacobian matrix. 

The performances and efficiency of the new method have been demonstrated through 
one dimensional shock tube problems as well as multidimensional problems. In terms of 
accuracy, the proposed method of solution greatly improves over the traditional single- 
state HLL approach in that it provides adequate and superior resolution of rotational 
discontinuities and a reduced numerical diffusion. Furthermore, it has been found that 
the new "HLLD" solver yields accuracies comparable to or even better than the Roe 
solver by retaining simplicity and greater ease of implementation over the latter. In terms 
of efficiency, the proposed HLLD solver is only ^ 20 — 25% slower than the HLL method 
but considerably faster (up to ~ 60% in some cases) than the Roe scheme. 

The author whishes to thank the support and hospitality received at the center of magnetic 
self-organization at the University of Chicago. 
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